BIO 356 Lab: Predator-Prey Dynamics
Rafael D’Andrea - Stony Brook University - Spring 2023
In this lab you will learn
Predators and prey regulate each other’s populations.
Predator-prey systems tend to be cyclic: predator and prey populations oscillate over time in response to each other.
This cycling can be stabilized by prey self-regulation, predator satiation, and prey switching.
A predator can promote coexistence of two competing prey species by limiting the growth of the superior competitor.
A keystone species is usually a predator with an outsize influence on the ecosystem. Their removal leads to a trophic cascade down the food web.
Release from predators may explain why some species become invasive when introduced to a new region.
library(tidyverse)
library(gridExtra) ## for plotting function grid.arrange()
library(deSolve) ## for predator-prey-resource ODE modeltheme_set(theme_bw())
theme_update(panel.grid = element_blank())
Ecology studies how organisms interact with the environment and other organisms. Of the different types of interactions between living organisms, extractive interactions are one of the most fundamental: in order to stay alive, grow, and reproduce, living organisms must constantly extract matter and energy from their surroundings. Plants and other photosynthetic and chemosynthetic organisms obtain their nutrients and energy directly from the abiotic environment, while all others rely on them for their own nutrition via trophic interactions. This simple fact of life gives rise to trophic chains (e.g. phytoplankton \(\rightarrow\) krill \(\rightarrow\) fish \(\rightarrow\) penguin \(\rightarrow\) seal \(\rightarrow\) orca) and food webs like the one below.
Figure: a food wweb in a temperate deciduous forest. Image from here.
Food webs are the means by which energy and matter are transferred across trophic levels, from primary producers (i.e. autotrophic organisms like plants and algae) to primary consumers to secondary consumers, etc. The biophysics of these transfers impose important constraints on life on Earth.
Q1: The British biologist Charles Elton, one of the founders of population and community ecology, noted in his influential book Animal Ecology (1927) that food chains rarely have more than five trophic levels. Based on content from lecture this week, what hypothetical explanation would you offer for this empirical observation?
Figure: Charles Elton (1900-1991)
The fact that the subsistence of an organism relies on the presence of another suggests important effects of trophic interactions on the dynamics of entire ecosystems. We would thus like to understand what kinds of effects these trophic interactions have on the populations involved. For that, we need a model. And in observance of the all-important KISS principle, we must start simple.
Basic predator-prey dynamics
Wolves hunting an elk at Yellowstone National Park. Photo by the National Park Service
Imagine a system consisting of a predator species and a prey species. The prey is the only source of food to the predator, and is not preyed on by any other predator. For example we can think of gray wolves and elk at Yellowstone National Park (this is a simplification as both species have trophic interactions with other park residents, although to a lesser extent than with each other). We will keep track of the two populations through annual censuses. We can write down a model describing how those populations are expected to change every year, as follows:
\[\begin{align} \tag{1} V_{t+1} = (b_V - d_V) V_t\\ \tag{2} P_{t+1} = (b_P - d_P) P_t \end{align}\]
Here, \(V\) denotes the prey (V for “victim”, so we can have different initials for predator and prey), and \(P\) the predator. This is the Nicholson-Bailey model of predator-prey dynamics, which is essentially a discrete-time version of the Lotka-Volterra equations. In its most general formulation, all of the rates above \(b_V, d_V, b_P, d_P\) depend on the prey and predator population.
However, let’s start simple. Let’s assume the prey have a constant per-capita rate of population growth \(b_V = b\) in the absence of the predator. That rate reflects the difference between the prey’s natural birth and death rates. This means that without the predator, each prey would always contribute \(b\) new individuals to next year’s prey population.
Q2. Without predators, what would happen to the prey population?
Predation causes additional per-capita mortality \(d_V\) on the prey. Let’s assume every year each predator kills a fixed proportion of the prey population. The number of prey each predator kills per year is then \(aV_t\). The constant of proportionality \(a\) is called the attack rate, and reflects the predator’s ability to find and catch prey. Therefore the total number of killings caused by this year’s predator population \(P_t\) will be \(aV_tP_t\). Comparing this to Equation (1), we conclude that \(d_V = a P_t\). Notice that unlike the natural growth rate \(b\) of the prey, this death-by-predation rate varies over time, as it depends on the size of the predator population.
Let’s also assume the predators have their own natural death rate per-capita death rate \(d_P = d\), which is a constant. This means that every year a constant fraction \(d\) of the predator population dies. Their birth rate \(b_P\), on the other hand, depends on their ability to convert captured prey into offspring. We can then write \(b_P = eaV_t\), where \(aV_t\) is the number of killings each predator makes per year and \(e\) is the conversion efficiency. The higher this efficiency, the fewer prey a predator must consume to generate one new offspring.
Putting this all together, our model then reads
\[\begin{align} \tag{3} V_{t+1} = (b - aP_t) V_t\\ \tag{4} P_{t+1} = (eaV_t - d) P_t \end{align}\]
Q3. Without prey, what would happen to the predator population?
Here’s the table of symbols in our model and their respective meanings:
| Symbol | Meaning |
|---|---|
| \(P_t\) | predator population size at year \(t\) |
| \(V_t\) | prey population size at year \(t\) |
| \(a\) | attack rate of the predator |
| \(b\) | natural per-capita rate of growth of the prey (also called intrinsic net growth rate) |
| \(d\) | natural per-capita death rate of the predator |
| \(e\) | conversion efficiency |
Here is the R code for our predator (P) - prey (V) model:
PV_Model =
function(
V0,
P0,
prey_intrinsic_net_growth_rate,
prey_carrying_capacity = 1e6,
predator_intrinsic_death_rate,
attack_rate,
conversion_efficiency,
final_year,
handling_time = 0
){
V_record = V = V0
P_record = P = P0
for(time in 1:final_year){
bv = prey_intrinsic_net_growth_rate * (1 - V / prey_carrying_capacity)
dv = attack_rate / (1 + attack_rate * handling_time * V) * P
bp = conversion_efficiency * attack_rate / (1 + attack_rate * handling_time * V) * V
dp = predator_intrinsic_death_rate
V = pmax(0, round((bv - dv) * V))
P = max(0, round((bp - dp) * P))
V_record = c(V_record, V)
P_record = c(P_record, P)
if(any(V == 0, P == 0)) break
}
return(
list(
parameters =
list(
prey_intrinsic_net_growth_rate = prey_intrinsic_net_growth_rate,
prey_carrying_capacity = prey_carrying_capacity,
attack_rate = attack_rate,
conversion_efficiency = conversion_efficiency,
predator_intrinsic_death_rate = predator_intrinsic_death_rate,
final_year = final_year
),
initial_conditions =
list(
V0 = V0,
P0 = P0
),
state =
tibble(
time = 0:time,
V = V_record,
P = P_record
)
)
)
}And here is the code to plot the PV model outcomes:
PV_Plot =
function(model){
plot_phase =
model$state |>
ggplot(aes(V, P, color = time)) +
geom_path() +
geom_point() +
theme(aspect.ratio = 1) +
scale_color_gradient(
low = 'royalblue1',
high = 'tan2'
)
plot_time =
model$state |>
pivot_longer(
-time,
names_to = 'species',
values_to = 'abundance'
) |>
ggplot(
aes(
time,
abundance,
group = species,
color = species
)
) +
geom_point() +
geom_line() +
theme(aspect.ratio = 1)
grid.arrange(
plot_time,
plot_phase,
nrow = 1
)
}We can calculate the equilibrium abundances of predator and prey, \(P^*\) and \(V^*\), by finding the values of \(P\) and \(V\) in Equations (3) and (4) such that \(V_{t+1}=V_t\) and \(P_{t+1}=P_t\). A little reflection should convince you that this will happen if the terms in parentheses equal 1, that is, if \(b - aP_t = 1\) and \(eaV_t - d = 1\). After some algebra, this works out to be \(P^* = \frac{b \, - \, 1}{a}\) and \(V^* = \frac{1 \, + \, d}{ea}\).
Q4. Suppose the prey have natural net growth rate \(b=2\), and the predators have natural death
rate \(d=0.1\). Suppose the attack rate
is 0.05, and the conversion efficiency is 0.1. Calculate the equilibrium
abundances \(P^*, V^*\). Calculate the
equilibrium abundances \(P^*\) and
\(V^*\), then set the initial
abundances \(P_0, V_0\) to those
values. Run the model and verify that these abundances do not change
over time. You can use the code below, replacing the ??
with the appropriate values. Show your plots.
PV_Plot(
model =
PV_Model(
V0 = ??,
P0 = ??,
prey_intrinsic_net_growth_rate = ??,
predator_intrinsic_death_rate = ??,
attack_rate = ??,
conversion_efficiency = ??,
final_year = 100
)
)The exercise above tells us that both the predator and prey populations have an equilibrium value, such that if they ever reach these populations at the same time, they will remain there forever. So this sounds like a sustainable system. But now let’s see if we should expect these populations to remain in equilibrium.
Q5. Rerun the model under the same parameters as above, but this time start the prey population at \(V_0 = 220\) and the predator population at \(P_0 = 21\). What happens? Feel free to check out what happens under several other initial conditions. What does this result say about the prospects of coexistence between the predator and the prey in real life, where random events may constantly cause small perturbations to their equilibrium abundances?
As we can see, this model is not optimistic about the prospects of sustainable predator and prey populations. However, there are many examples in nature of long-term coexistence between predator and prey. This suggests there is one or more important elements missing from our model that allows this coexistence in nature.
Note: Of course, models are always simplifications of reality, and therefore are always missing elements. Indeed, one of the principles of science and statistics is that All models are wrong, but some are useful (George Box). In other words, the question is not whether the model is right but whether it is informative. From the finding that one cannot get a sustainable predator-prey system under the assumptions of our model, we learn that some of the missing elements are critical for predator-prey coexistence. That is what we turn to next.
Prey self-regulation
As we’ve seen in lecture, one critical element to a sustainable population is self-regulation. This element is missing from the prey dynamics: without predators, the prey grow unchecked to infinity. This is obviously unrealistic, and a good candidate for the important missing element. So let’s introduce it to our model.
Recall from earlier in the course that the simplest form of self-regulation is a quadratic term in the equation for population growth. We saw this when we studied logistic growth, \(\frac{dN}{dt} = rN\left(1 - \frac{N}{K}\right)\). The quadratic term is the second term in the parenthesis. The quantity \(K\) in that term is the carrying capacity, which is the equilibrium value of the population in that model. Let us then add self-regulation to the prey via the introduction of a carrying capacity. Looking at the code, you will notice that the carrying capacity is already present in the model, but its default value is set to 1 million prey, which is so large as to be indistinguishable from infinity for practical purposes.
Now let’s set the carrying capacity to a reasonable value, say, 1000 individuals, and see what happens:
PV_Plot(
model =
PV_Model(
V0 = 220,
P0 = 21,
prey_intrinsic_net_growth_rate = 2,
prey_carrying_capacity = 1000,
predator_intrinsic_death_rate = .1,
attack_rate = .05,
conversion_efficiency = .1,
final_year = 100
)
)
We still depart from the equilibrium values \(P^*, V^*\), never to return. As before, \(P^*, V^*\) is an unstable equilibrium. However, now the predator and prey species can sustain long-term populations! Both populations are bound between minimum and maximum values, and oscillate between them indefinitely.
How did the introduction of a finite carrying capacity to the prey cause such an impact?
As you saw in question Q4, when the prey population decreased due to predation from the initial extra predator, the predator population dropped to low numbers. In response, the prey population skyrocketed. This ultimately led the predator population to catch up with the prey and skyrocket as well, to more than twice its initial value. The resulting extra predation turned the balance between growth and mortality of the prey, and then its population crashed before it could recover. The fundamental problem here is that there is a cycle of feedback escalation between the predator and the prey.
On the other hand, the introduction of prey self-regulation reflected in its carrying capacity tamps down the rise in the prey population following the initial decrease in predation. This braking stops the escalation, and the predator population never reaches dangerously high numbers for the prey. Eventually, both populations become locked in a never-ending loop called a limit cycle.
Importantly, the cycle will be eventually reached even under different initial conditions, as you can verify by rerunning the model with different initial populations. This is critical because it means that random events such as catastrophic weather would have only temporary effect, as the populations would eventually fall back on their limit cycle.
Q6. What happens if we increase the carrying capacity to 2000? Which phenomenon discussed in class does this outcome exemplify? Given what you’ve seen in the Lab so far, how would you explain the occurrence of the phenomenon?
Predator satiation
We have been assuming that every year each predator eats a fixed proportion of the prey population. This means there is no limit to the number of prey individuals a single predator can eat in a fixed period of time. More realistically, predators will get satiated after eating a certain amount of food, and their bodies cannot process more food until they’ve had time to digest it. The minimum time required between predation events is commonly called handling time, and includes the time to find, catch, ingest, and digest prey. Because of handling time, each predator will have a maximum number of prey it can eat in a year even under an infinite prey population. From the prey’s perspective, the per capita risk of predation will be inversely proportional to the prey population.
We can introduce these considerations to our model by writing the number of prey that each predator eats in a year as \(\frac{aV}{1 + ahV}\), where \(h\) is the handling time. Notice that if \(h = 0\) we recover the original consumption amount \(aV\), which is how many prey the predators would eat under no time constraints. Our code in fact already includes the handling time, but we have been using the default value of \(h = 0\).
Note: Some algebra reveals that under the above formulation, the maximum number of prey that a predator can eat in a year is \(n_{\max} = 1/h\), but they will only be able to catch this many prey if the prey population is very large (for the mathematically inclined, the condition is \(V \gg \frac{1}{ah}\)). Under those circumstances, the annual risk to an individual prey of being caught by a predator is \(n_{\max} / V\).
Now let’s set the handling time to \(h = 0.02\) (meaning that a predator can eat at most 50 prey per year), and see what it does to our stabilized system above.
PV_Plot(
model =
PV_Model(
V0 = 220,
P0 = 21,
prey_intrinsic_net_growth_rate = 2,
prey_carrying_capacity = 1000,
predator_intrinsic_death_rate = .1,
attack_rate = .05,
conversion_efficiency = .1,
final_year = 100,
handling_time = 0.02
)
)
We still have a limit cycle. However, this is a tighter cycle: both predator and prey populations oscillate within a narrower range than before. This suggests that predator satiation is further stabilizing our system, presumably by further dampening the feedback escalation between prey and predator.
Q7. Set the code back to the settings from Q5 and try raising the handling time to \(h = 0.03\). What happens? How do you interpret the contribution of predator satiation to the prospects of a stable equilibrium between predators and prey?
Prey switching
What if the predator had an alternative food source?
Suppose we introduce another prey species to our park. The second prey has a similar size and nutritional value to the predator, which therefore has a similar attack rate and conversion efficiency as for the first prey species. However, the two prey species reside in different habitats in the park, and the predator spends most of its time stalking whichever prey is currently easier to catch. As a result, each year the predator only attacks whichever prey species is more abundant. This is a simplified version of prey switching behavior, on which there is extensive modeling literature (for one such study on the predatory behavior of wolves on elk and bison in Yellowstone National park, check out this article) and a number of observational studies and experiments (examples here and here).
One could argue that the introduction of the second prey could further stabilize the system because the predator always has food in quantity while overhunted prey is allowed to recover. Let’s test this idea on our model. The code for the 1-predator 2-prey version of the model is provided below, along with code to plot results.
PVV_Model =
function(
V0,
P0,
prey_intrinsic_net_growth_rate,
prey_carrying_capacity = c(1e6, 1e6),
predator_intrinsic_death_rate,
attack_rate,
conversion_efficiency,
final_year,
handling_time = 0,
predation_mode = 'full_switch'
){
V_record = V = V0
P_record = P = P0
for(time in 1:final_year){
if(predation_mode == 'full_switch') alpha = attack_rate * (V == max(V))
if(predation_mode == 'proportional_switch') alpha = attack_rate * V / sum(V)
bv = prey_intrinsic_net_growth_rate * (1 - V / prey_carrying_capacity)
dv = alpha / (1 + alpha * handling_time * V) * P
bp = sum(conversion_efficiency * alpha / (1 + alpha * handling_time * V) * V)
dp = predator_intrinsic_death_rate
V = pmax(0, round((bv - dv) * V))
P = max(0, round((bp - dp) * P))
V_record = rbind(V_record, V)
P_record = c(P_record, P)
if(any(V == 0, P == 0)) break
}
return(
list(
parameters =
list(
prey_intrinsic_net_growth_rate = prey_intrinsic_net_growth_rate,
prey_carrying_capacity = prey_carrying_capacity,
attack_rate = attack_rate,
conversion_efficiency = conversion_efficiency,
predator_intrinsic_death_rate = predator_intrinsic_death_rate,
final_year = final_year
),
initial_conditions =
list(
V0 = V0,
P0 = P0
),
state =
tibble(
time = 0:time,
V1 = V_record[, 1],
V2 = V_record[, 2],
P = P_record
)
)
)
}Code for plotting:
PVV_Plot =
function(model, plots = 'time_series'){
plot_phase_VV =
model$state |>
ggplot(aes(V1, V2, color = time)) +
geom_path() +
geom_point() +
theme(aspect.ratio = 1) +
scale_color_gradient(
low = 'royalblue1',
high = 'tan2'
)
plot_phase_PV1 =
model$state |>
ggplot(aes(V1, P, color = time)) +
geom_path() +
geom_point() +
theme(aspect.ratio = 1) +
scale_color_gradient(
low = 'royalblue1',
high = 'tan2'
) +
scale_color_gradient(
low = 'royalblue1',
high = 'tan2'
)
plot_phase_PV2 =
model$state |>
ggplot(aes(V2, P, color = time)) +
geom_path() +
geom_point() +
theme(aspect.ratio = 1) +
scale_color_gradient(
low = 'royalblue1',
high = 'tan2'
)
plot_time =
model$state |>
pivot_longer(-time, names_to = 'species', values_to = 'abundance') |>
ggplot(aes(time, abundance, group = species, color = species)) +
geom_line(linewidth = 1) +
theme(aspect.ratio = 1)
if(plots == 'all'){
grid.arrange(
plot_time,
plot_phase_VV,
plot_phase_PV1,
plot_phase_PV2,
nrow = 2
)
}else plot_time
}We can now run the model. Let’s assume the second prey species \(V_2\) has intrinsic net growth rate 2.1 and carrying capacity 750, compared to 2 and 1,000 respectively for prey species \(V_1\). Let’s say we introduce 10 individuals of prey species 2. Notice that now in the model call, the prey parameters are vectors, with the first element referring to the original prey species \(V_1\) and the second element referring to the introduced species \(V_2\).
model =
PVV_Model(
V0 = c(220, 10),
P0 = 21,
prey_intrinsic_net_growth_rate = c(2, 2.1),
prey_carrying_capacity = c(1000, 750),
predator_intrinsic_death_rate = .1,
attack_rate = .05,
conversion_efficiency = .1,
final_year = 50,
handling_time = 0.02,
predation_mode = 'full_switch'
)
PVV_Plot(
model = model,
plots = 'all'
)
Interestingly, the addition of the second prey destabilized the system! All three species coexist in the park, but now the population oscillations are higher than before under similar parameters, particularly for the prey species! Also, the oscillations are now less predictable, as we can see by examining the phase plots between the predator and each prey, and between the two prey.
Q8. Rerun the 1-predator-2 prey model with predator handling
time 0.03 and 0.04. Show your plots. Check the table of population sizes
over time for the ten last years of the simulation by running the line
model$state |> tail(10). What happens to the population
of the predator and the two prey as you raise the handling time from
0.02 to 0.03 to 0.04?
Let’s insist on this idea a bit further. Perhaps rather than an all-or-nothing full switch between prey species from one year to the next, the predator may actually attack prey in proportion to how abundant they are. For example, if prey species 1 is currently 60% of both prey populations and prey 2 makes up the other 40%, then the predator will invest 60% of its preying activity on prey 1 and the remaining 40% on prey 2.
Q9. Rerun the model with the same prey and predator
parameters as above (use handling time 0.04), this time setting the
argument predation_mode from 'full_switch' to
'proportional_switch'. What happens now?
Under this second prey-switching strategy, the predator has a much more refined response to prey abundance, and is thus able to exert a more fine-tuned density-dependent control over the population of the prey. Consequently, the dynamics are further stabilized.
Those two prey-switching strategies we just examined represent the ends of a spectrum. In real life, prey-switching behavior will likely fall somewhere in between. However, as we have just seen, the details of the predator’s behavior can make all the difference. Overall, predator-prey models are notorious for stability issues, and a great deal of ecological research has been dedicated to finding sources of stabilization.
Trophic cascades
Finally, consider a scenario with not two but three trophic levels: a predator, two prey species, and a resource. Both prey species consume the resource, and both are eaten by the predator. For concreteness, consider a pond populated by two species of phytoplankton – a diatom and a green alga – and a copepod (a small crustacean) which feeds on both the diatom and the alga. The phytoplankton in turn grow on nitrate, which flows into the pond via a tributary stream.
Figure: Stigeoclonium, a green alga, by Kristian Peters. Diatoms, by
Prof. Gordon T. Taylor, Stony Brook University. A copepod, by Otto
Larink - Plankton Net.
Given the microscopic nature of these organisms and the abiotic resource, it makes more sense to model their biomass than their head count. In addition, the dynamics of nitrate flow and phytoplankton/zooplankton growth are better described as a ongoing process than an annual process. As such, we will now use a continuous-time model. One of the main differences from the discrete-time scenario is that now our species respond immediately to changes in each other’s abundances. In practice, these systems are more easily stabilized.
Let’s write down our model:
\[ \frac{dR}{dt} = s - c_1 \, R \, V_1 - c_2 \, R \, V_2 \\ \frac{dV_1}{dt} = (e_{v_1} \, c_1 \, R - a_1 \, P-d_{v_1}) \, V_1 \\ \frac{dV_2}{dt} = (e_{v_2} \, c_2 \, R - a_2 \, P-d_{v_2}) \, V_2 \\ \frac{dP}{dt} = (e_p \, a_1 \, V_1 + e_p \, a_2 \, V_2 - d_p) \, P \]
In words, the resource \(R\) flows in at a constant rate \(s\) and is depleted via consumption by both phytoplankton species \(V_1\) (the diatom) and \(V_2\) (the alga), according to their consumption rates \(c_1\) and \(c_2\), reflecting each species’ ability to consume and process nitrate. Both the diatom and the alga grow as a result of consuming nitrate, and their biomass is depleted both via natural death at rates \(d_{v_1}\) and \(d_{v_2}\), and via predation by the copepod. The copepod eats the diatom at rate \(a_1\) and the alga at rate \(a_2\), reflecting its preference for each prey, and has its own intrinsic mortality \(d_p\). The efficiencies of the prey \(e_{v_1}, e_{v_2}\) and predator \(e_p\) determine the rate at which they grow as they eat.
The symbols in the model and their meanings are laid out in the table below
| Symbol | Meaning |
|---|---|
| \(R\) | concentration of the resource – nitrate |
| \(V_1\) | population size of prey species 1 – diatom |
| \(V_2\) | population size of prey species 2 – alga |
| \(P\) | population size of the predator – copepod |
| \(s\) | supply rate of nitrate |
| \(c_{1(2)}\) | consumption rate of nitrate by prey 1 (2) |
| \(e_{v_{1(2)}}\) | conversion efficiency of prey 1 (2) |
| \(e_p\) | conversion efficiency of the predator |
| \(a_{1 (2)}\) | attack rate of the predator on prey 1 (2) |
| \(d_{v_{1(2)}}\) | mortality of prey 1 (2) |
| \(d_p\) | mortality of the predator |
The following R code encapsulates this model.
PVVR_Model =
function(
R0,
V0,
P0,
resource_supply,
consumption_rates,
predation_rates,
prey_mortality,
predator_mortality,
prey_efficiency,
predator_efficiency,
final_year,
time_step
){
PVVR = function(t, state, parameters){
with(
as.list(c(state, parameters)), {
dRdt = s - c1 * V1 *R - c2 * V2 *R
dV1dt = (ev1 * c1 * R - a1 * P - dv1) * V1
dV2dt = (ev2 * c2 * R - a2 * P - dv2) * V2
dPdt = (ep * a1 * V1 + ep * a2 * V2 - dp) * P
list(c(dRdt, dV1dt, dV2dt, dPdt))
}
)
}
parameters =
c(
s = resource_supply,
c = consumption_rates,
a = predation_rates,
dv = prey_mortality,
dp = predator_mortality,
ev = prey_efficiency,
ep = predator_efficiency
)
state = c(R = R0, V1 = V0[1], V2 = V0[2], P = P0)
times = seq(0, final_year, by = time_step)
out = ode(y = state, times = times, func = PVVR, parms = parameters)
return(
list(
parameters =
list(
resource_supply,
consumption_rates,
predation_rates,
prey_mortality,
predator_mortality,
prey_efficiency,
predator_efficiency,
final_year = final_year
),
initial_conditions =
list(
R0 = R0,
V0 = V0,
P0 = P0
),
state = out
)
)
}Because this model has four “dimensions” (the predator, two prey, and the resource), it contains more parameters. We will be focusing on the parameters that distinguish the diatom from the green alga, namely their respective affinities for the resource, the preference of the copepod for eating each of them, and their intrinsic mortalities.
Let’s consider the following parameters for the two phytoplankton species:
\[c_1 = 0.07 \hspace{1cm} c_2 = 0.03 \hspace{1cm} a_1 = 0.008 \hspace{1cm} a_2 = 0.002 \hspace{1cm} d_{v_1} = 0.02 \hspace{1cm} d_{v_2} = 0.13\]
Because of the continuous-time nature of the model, the numerical values of the parameters are less easy to interpret than before, but this does not concern us as we are interested in the comparison between the two species.
Q10. Based on the parameter values above, which phytoplankton species a) Is the faster consumer of nitrate? b) Is the copepod’s preferred prey? c) Has the lowest mortality?
After setting other parameters and initial abundances, now measured as concentration in the water, we are ready to run the model and see what happens:
model =
PVVR_Model(
R0 = 1000,
V0 = c(100, 40),
P0 = 20,
resource_supply = 1e4,
consumption_rates = c(.07, .03),
predation_rates = c(.008, .002),
prey_mortality = c(.02, .13),
predator_mortality = .01,
prey_efficiency = c(.01, .01),
predator_efficiency = .01,
final_year = 1500,
time_step = 1
)
plot =
model$state[, 1:5] |>
as_tibble() |>
rename(predator = 'P', resource = 'R', `prey 1` = V1, `prey 2` = V2) |>
pivot_longer(-time, names_to = 'species', values_to = 'abundance') |>
ggplot(aes(time, abundance)) +
geom_line(color = rgb(153 / 255, 0, 0)) +
facet_wrap(~species, scales = 'free') +
theme(aspect.ratio = 1)
plot
After wide initial oscillations, the concentration of the resource and the abundance of all species relaxed to a positive equilibrium.
Looking closely, we notice that at first prey 1 spikes in abundance and depletes the resource, leading to the almost extinction of prey 2. However, the predator responds to such a surplus of its favorite prey and its abundance quickly rises. As a result, the abundance of prey 1 quickly falls back down, allowing the resource to go back up, and ultimately also prey 2. In the end, all four reach a steady concentration.
An equilibrium consisting of fixed abundances, also called a point equilibrium (to be contrasted with cyclic equilibria like the ones we’ve seen in this lab), is stable if the system converges towards it from its vicinity. In other words, if the current abundances of the resource, prey, and predator are not too far from their final equilibrium abundances, they eventually will reach these abundances. In this case the point equilibrium is called a point attractor. Its opposite, a point repellor, is an equilibrium where abundances move away from it if they’re not exactly at it. You saw an example of a point repellor in Q4.
Q11. Check that the equilibrium above is a point attractor by rerunning the code under different initial abundances, and verifying that the same final abundances are reached. Show your code and one of your plots.
Q12. What happens if we removed the predator from the ecosystem? Try to predict the outcome first, then run the model with initial predator abundance set to 0. Explain the final abundances of the resource and each prey in terms of the parameters of the prey.
The outcome of Q12 touches on several important topics in community ecology and conservation biology: trophic cascades, keystone species, and the enemy release hypothesis. It also touches on another topic of great importance to community ecology: coexistence theory. Notice that the prey can coexist when both the resource and predator are present, but not in the absence of the predator (or the resource for that matter). Prey 1 is better at procuring the resource, but Prey 2 is better at avoiding the predator. This tradeoff in how each species handles the two limiting factors in our model is the key to their coexistence. We will turn our attention to the very important topic of species coexistence in our next module, Competition.
Q13. Noting what happened to the abundance of the resource upon the removal of the predator, explain the concept of a trophic cascade.
Keystone species
The removal of the copepod and subsequent extinction of one of the prey and depletion of the resource is reminiscent of Robert Paine’s famous 1960s experiments in the intertidal zone at Neah Bay in Washington State. Paine observed that the experimental removal of the purple sea star (Pisaster ochraceus) led to drastic alteration of the ecosystem, as one of the sea star’s prey, the California mussel (Mytilus californianus), ultimately dominated the previously diverse habitat. A similar example was a natural experiment following intensive hunting of sea otters for furs in the north Pacific coast in the 19th century: with otter numbers low or nonexistent at places, sea urchins were left without their main predator. Their population exploded, resulting in overgrazing and subsequent decimation of kelp forests, along with other species dependent on kelp (further details can be found here).
The purple sea star and the sea otter are considered keystone species, defined as species with an outsize impact on their ecosystem. As such, they are usually predators whose removal leads to a trophic cascade with consequences further down the food web.
Figure: Two keystone species. Purple sea star, by D. Gordon E. Robertson. Sea otter, by Marshal Hedin
Q14. Describe another example of keystone species. If you use a source, make sure to cite it.
Q15. A key characteristic of keystone species is that their influence on the ecosystem is disproportional to their abundance or biomass relative to other species. With this in mind, explain why keystone species are typically apex predators. Hint: Think about energy flow between trophic levels.
The enemy release hypothesis
When a non-native species is introduced to a new habitat, it may become invasive. An invasive species spreads upon introduction and becomes abundant to the point of altering their new environment, often adversely. Invasive plants typically grow faster, disperse faster, or otherwise outcompete native plants in some regard.
Ecologists have identified many potential reasons why some species become invasive. Chief among them is the enemy release hypothesis, which posits that in the new environment the species lacks the predators and parasites that kept its population in check in its native range.
Many empirical studies corroborate this hypothesis. For example, Wolfe (2002) found that native European populations of the white campion, Silene latifolia, are 17 times more likely to be damaged by aphids, fungal diseases, flower herbivory, and fruit predation, than their introduced North American couterparts. The author highlights that two of the plant’s main enemies, the seed predator insect Hadena rivularis and the anther smut fungus Microbotryum violaceum, are extremely common in the native range of the plant but absent in N. America.
Q16. Explain how the outcome of Q12 relates to the enemy release hypothesis.
Q17. What kind of biological control of pests and invasive species is logically suggested by the enemy release hypothesis and the results in Q12?